rm(list = ls())
setwd('E:\\17_HuaDong\\teach\\MyLecture\\R\\code')
library(foreign)
library(magrittr)
library(plm)
nlswork <- read.dta('DID.dta')
summary(nlswork)
# 生成控制变量的平方项
nlswork$age2 <- nlswork$age^2
nlswork$ttl_exp2 <- nlswork$ttl_exp^2
nlswork$tenure2 <- nlswork$tenure^2

# 生成时间和政策虚拟变量
nlswork$dt <- 0
nlswork$dt[nlswork$year >= 77] <- 1 # 1977年以后才实施政策
nlswork$treated <- 0
nlswork$treated[nlswork$idcode > 2000] <- 1 # 大于2000的id才被处理

# 只有两期，加入dt,treated，如下回归即可
lm(ln_wage~dt+treated+dt:treated+grade+age+age2+ ttl_exp+ ttl_exp2+
     tenure+ tenure2+ not_smsa+ south+ race,data = nlswork) %>% summary()
# 多期的话，要加入时间和个体固定效应，如下回归
plm(ln_wage~dt:treated+grade+age+age2+ ttl_exp+ ttl_exp2+
     tenure+ tenure2+ not_smsa+ south+ race,data = nlswork,
    index = c('idcode','year'),effect = 'twoways',model = 'within') %>% summary()

# 稳健性分析
## 共同趋势检验：以政策实施前的样本，随机寻找一个伪实施年
nlswork$dt72 <- 0
nlswork$dt72[nlswork$year >= 72] <- 1 # 假设1972年以后才实施政策
plm(ln_wage~dt72:treated+grade+age+age2+ ttl_exp+ ttl_exp2+
      tenure+ tenure2+ not_smsa+ south+ race,data = nlswork[nlswork$year<77,],
    index = c('idcode','year'),effect = 'twoways',model = 'within') %>% summary()
## 控制组不受政策影响：对照两个控制进行前述分析
nlswork$PseTreated <- 0
nlswork$PseTreated[nlswork$idcode > 1000 & nlswork$idcode < 1600] <- 1 # 大于2000的id才被处理
plm(ln_wage~dt:PseTreated+grade+age+age2+ ttl_exp+ ttl_exp2+
      tenure+ tenure2+ not_smsa+ south+ race,data = nlswork[nlswork$idcode < 2000,],
    index = c('idcode','year'),effect = 'twoways',model = 'within') %>% summary()


# 动态DID
## 生成D_+1等超前、滞后变量
for (i in unique(nlswork$year)) {
  MyStr <- paste('nlswork$year',as.character(i),'<- 0',sep = '')
  eval(parse(text = MyStr))
  MyStr <- paste('nlswork$year',as.character(i),'[nlswork$year==',as.character(i),'] <- 1',sep = '')
  eval(parse(text = MyStr))
  MyStr <- paste('nlswork$treated',as.character(i),'<- nlswork$year',
                 as.character(i),'*nlswork$treated',sep = '')
  eval(parse(text = MyStr))
}
## 1968-1975年，政策实施前的不应显著:本例说明预期效应存在
MyStr <- paste('treated',as.character(sort(unique(nlswork$year))[1:7]),sep = '') %>%
  paste(.,collapse = '+') %>% 
  paste('fml <-ln_wage~',.,'+grade+age+age2+ ttl_exp+ ttl_exp2+tenure+ tenure2+ not_smsa+ south+ race',sep ='')
eval(parse(text = MyStr))
plm(fml,data = nlswork,
    index = c('idcode','year'),effect = 'twoways',model = 'within') %>% summary()
## 滞后影响
MyStr <- paste('treated',as.character(sort(unique(nlswork$year))[7:14]),sep = '') %>%
  paste(.,collapse = '+') %>% 
  paste('fml <-ln_wage~',.,'+grade+age+age2+ ttl_exp+ ttl_exp2+tenure+ tenure2+ not_smsa+ south+ race',sep ='')
eval(parse(text = MyStr))
plm(fml,data = nlswork,
    index = c('idcode','year'),effect = 'twoways',model = 'within') %>% summary()
